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Abstract 

Sufficiently thin films of type-I superconductor in a perpendicular magnetic field exhibit a tri- 
angular vortex lattice, while thick films develop an intermediate state. To elucidate what happens 
between these two regimes, precise numerical calculations have been made within Ginzburg-Landau 
theory at n = 0.5 and 0.25 for a variety of vortex lattice structures with one flux quantum per unit 
cell. The phase diagram in the space of mean induction and film thickness includes a narrow wedge 
in which a square lattice is stable, surrounded by the domain of stability of the triangular lattice 
at thinner films/lower fields and, on the other side, rectangular lattices with continuously varying 
aspect ratio. The vortex lattice has an anomalously small shear modulus within and close to the 
square lattice phase. 

PACS numbers: 74.25.Uv, 74.20.De, 74.78.-w 
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I. INTRODUCTION 



Thin films of a bulk type-I superconductor subject to a perpendicular magnetic field can 
behave like bulk type-II superconductors, in that they develop a vortex lattice. The pioneer- 
ing theoretical treatments in the area, by Tinkhami and Maki 2 applied Ginzburg-Landau 
(GL) theory in the vicinity of the critical field where the order parameter vanishes. Pearl's^ 
treatment of isolated vortices within London theory, appropriate in the low-field limit, shows 
that vortices in a sufficiently thin film have a long-range repulsion. This repulsion should 
lead to the development of a triangular vortex lattice at small fields, and it was shown 
by Lasher 5 that a triangular vortex lattice is also favored near the upper critical field for 
sufficiently thin films. 

On the other hand, sufficiently thick films of type-I superconductors exhibit the interme- 
diate state. How does the magnetic flux structure evolve from triangular vortex lattice to 
intermediate state with increasing film thickness? In other words, what is the equilibrium 
phase diagram for magnetic flux structures, as a function of film thickness and magnetic 
field? In real superconductors the picture will necessarily be complicated by disorder and 
anisotropy, but the question is interesting enough within GL theory. Lasher- was the first 
to address it. He showed that, within linearized GL theory, between the triangular vortex 
lattice and the intermediate state there were a large number of distinct vortex phases, with 
the triangular lattice first replaced by a square vortex lattice. Some years later Callaway^ 
pointed out that Lasher had not considered the most general Abrikosov-type solutions to 
the linearized first GL equation, and he carried out a comprehensive analysis of periodic 
vortex arrays, valid close to the upper critical field. 

While theoretical analyses of the magnetic flux structure phase diagram have been re- 
stricted to the highest possible fields, interesting experimental results have appeared at very 
low fields. Hasegawa et air- applied electron holography to examine the magnetic field in 
the space above flux structures in Pb films. They found evidence for vortices with more 
than one flux quantum (which they denoted MQF-A) as well as flux structures that seemed 
more likely to be associated with normal regions of finite cross-section (which they denoted 
MQF-B). "Multiply-quantized" (also known as "giant") vortices are known to arise in various 
circumstances. Holes in a superconductor parallel to the field trap vortices with greater 
fluxoid number as their radii increase;- the repulsion of vortices from a film edge can lead to 
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the formation of a equilibrium giant vortex in the center of a small, thin disk;- metastable 
giant vortices develop in field-cooling of small cylinders. 10 None of these seem relevant to 
the experiment of Hasegawa et ai, and the search for stable lattices of multiply-quantized 
vortex in the phase diagram for type-I films was one of the motivations for the present work. 

The magnetic flux structure phase diagram is also of fundamental interest, because vor- 
tices in a film with GL parameter k < 1/ v^2 comprise an archetypical system with competing 
interactions. It is not obvious whether the vortex structures found by Lasher and Callaway 
at intermediate thicknesses survive on reducing the magnetic field, and attempts to experi- 
mentally observe such structures would benefit from theoretical guidance. 

It is also noteworthy that a bulk GL superconductor with k = 1/ \[2 and at the critical 
field exhibits massive (in fact complete) degeneracy with respect to vortex configurations.— 
Luk'yanchuk 12 has carried out a thorough analysis of corrections to the GL functional, 
together with deviations of k and the magnetic field from their critical values, in breaking 
the degeneracy. He noted that demagnetization effects also break the degeneracy, but did no 
calculations along those lines. A film geometry corresponds to maximum demagnetization, 
so it may be interesting to compare the vortex phase diagram for films with /c«l/ \/2 with 
the phase diagrams that follow from the analysis by Luk'yanchuk. 

In this paper we take the first steps towards filling out the magnetic flux structure phase 
diagram for the minimal model, isotropic Ginzburg-Landau theory, of thin film type-I super- 
conductors. The competition between various phases is delicate at the upper critical field, 
so accurate free-energy calculations for different flux structures are necessary. Consequently, 
we have followed the approach pioneered by Brandt for vortex lattices in bulk 13 and, more 
recently, thin film^ GL superconductors. The squared magnitude of the order parameter, 
the supervelocity, and the magnetic field are represented as linear combinations of appro- 
priate basis functions. The GL equations then become a set of nonlinear equations for the 
coefficients, which are solved by iteration. 

Sec. [Til describes the computational method in more detail. Brandt's papers are quite 
explicit, so we may be brief, and highlight the modest changes we made to Brandt's algorithm 
for thin films superconductors. The present calculations are restricted to magnetic flux 
structures consisting of singly-quantized vortices in periodic structures with one vortex per 



unit cell. We have found that the functional form chosen for the magnetic field in Ref. 
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limits the accuracy of the magnetic field and consequently the free energy, and we offer a 
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correct and computationally convenient alternative. Sec. II III presents the principal results of 
the calculations, which are based on evaluating the free energy for a large number of points 
in the space of vortex lattices structures, film thickness, and applied magnetic fields (or, 
equivalently, mean inductions). Phase diagrams, free energy densities, and vortex lattice 
shear moduli are given for k = 0.5 and 0.25. Other values of the GL parameter could 
have been considered but the calculations become significantly more challenging at smaller 
values of k; and with results for just two values, some trends with variation of k may be 
deduced. Sec. [TV] offers various decompositions of the free energy density to facilitate the 
physical interpretation of the flux structure phase diagram. In Sec. |V] we summarize the 
results and note their limitations, indicate some directions for future theoretical work, and 
offer suggestions for experiments. 



II. COMPUTATIONAL METHOD 



In order to avoid unnecessary repetition of material presented in Ref. 14, we will start by 
presenting only as much of it needed to make our further developments intelligible. 

The standard reduced units are employed, in which lengths are in units of the penetration 
depth A, energy densities are in units of fioH 2 (with H c the thermodynamic critical field) 
and magnetic induction is in units of \/2^qH c . Note that in these units the upper critical 
mean induction is k. We consider infinite films with — d/2 < z < d/2. The GL free energy 
can be expressed in gauge-invariant form (rather than in terms of the GL order parameter 
if) = fe lip and vector potential A) by writing it in terms of the square of the order parameter, 
lo = f 2 , the supervelocity, 

Q = A-K~ 1 V<p (1) 
and the deviation from mean induction, 



B-zB 



(2) 



We let S denote both the unit cell area and the unit cell itself, depending on context; for 
the former, S = B/<f>o with $o — 2tt/k. The free energy per unit volume of superconductor 
referenced to the normal state is 



Sd 



dxdy 



d/2 



dz 



d/2 



1 2 I Vo;| ^ 2 



Sd 



dxdy I dz[B 2 

d/2 



B 2 } 
(3) 



where the contribution of the first two terms in the first integral is the condensation free 
energy F con d, that of the next two terms is the kinetic energy of the supercurrent Fki n , that 
of the last term is the internal field energy F mag , and that of the second integral is the stray 
field energy F stray . In order to determine the phase diagram we will compare the minimum 
F for different vortex lattice structures with the same value of B (and hence S). For a 
particular vortex lattice structure, minimizing the GL free energy with respect to variations 
in the order parameter yields the first GL equation 

2k 2 \ ' 2u 



The second GL equation is 

V x B = -cuQ (5) 
which is identical to Ampere's law in reduced units because the supercurrent j is 

j = -uQ (6) 

A key step in Brandt's approach is to decompose the supervelocity as 

Q = Qa + q (7) 

where is the supervelocity of the Abrikosov solution corresponding to the given vortex 
lattice, which satisfies 



V x Q A 



B - $ ^ (l\L - Rvortex)] Z (8) 



where = (x,y), S2 is the two-dimensional Dirac delta function, and the sum runs over all 
vortices. Inside the film, 

b = V x q (9) 

With these definitions and relations, the problem is to determine ou, q and b that minimize 
the free energy. Brandt's Ansatz for these fields is as follows: 

w(r)= Y a K ± K z [1 - cos K ± • rj cos (K z z) (10) 
k±,k z 

q(r)= Y b K±K z - sin Kj_ ■ r ± cos K z z (11) 
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b z (r) = &k ± a' z cos K ± ■ r ± cos K z z (12) 

bj_(r) = bjt ±Kz ^ * sin K ± • r ± sin if z 2: (13) 

Here is the set of reciprocal lattice vectors for the vortex lattice and K z = (2ir/d) n with 
n running over the whole numbers. Several features of this Ansatz are worth noting. Only 
two sets of expansion coefficients, &~k_^k z and &k ± k z; are required because b and q are linked 
by ([9]). The periodicity of the u combined with the quadratic behavior of u near the vortices 
suggests the form of expansion for the r_i_ dependence in ( fTUl) . while the boundary condition 
for the order parameter at a superconductor-insulator interface makes the cosine expansion 
natural for the z dependence. Eq. tTTTT) leads to supercurrents with, as one would anticipate, 
only in-plane components, as well as with the appropriate periodicity and behavior near 
vortex lines. The motivation for the z dependence of the expansions for b and q is that q 
and b z are even functions of z while is an odd function. 

Inserting f fTOl) and (llljl into the first GL equation and applying orthogonality of trigono- 
metric functions leads to coupled nonlinear equations for the expansion coefficients an ± K z 
and b-K ± K z which can be readily cast in the form of equations for the clk ± k z suitable for 
solution by iteration: see Eq. (T2~T|) below. More equations must come from the second 
GL equation inside the film, together with V x b = outside the film and the boundary 
conditions on the induction. The induction above the film satisfies 

B z = B + b k ± cos K ± ■ r ± e- K ^ z - d/2) (14) 



K 



B_l = $>K,|7 sinK ^ • r±e~ K ^- d/2) (15) 
and the continuity-of-i^ boundary condition may be expressed as 

b s K± = J2 b ^ z cosdK z /2 (16) 

It is convenient to derive the equations for the expansion coefficients by direct minimization 
of the free energy (including the stray field energy) with respect to the b^ ± K z , which leads 
to Eqs. (19) through (23) of Ref. [ijl which we will not reproduce here. 

In order to carry out a calculation of the expansion coefficients it is necessary to trun- 
cate the expansion, setting ci~k ± k z and b~K ± K z to zero for K±K Z outside some range. It is 
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also necessary to approximate the integrals that appear in the iteration equations as finite 
sums. Those integrals arise from applying orthogonality relations and, ideally, the coefficient 
truncation and numerical integration could be done consistently, so that the trigonometric 
functions retained in the expansion are orthogonal with respect to the numerical integration. 
This is done naturally for the z coordinates of the integration, by making the simplest choice 
of uniform spacing. In the xy plane Brandt employs a rectangular grid for integration but 
a circular domain for the allowed Ki values. Though a rectilinear domain for Ki would be 
more consistent we have followed Brandt's choice, on the grounds that when Kj_ is large the 
expansion coefficients ought to be small. 

What is there to object to in the method described above? In brief, Eq. (TTTj) (and its 
corollaries Eqs. ( lT2j) and f)13p ) impose periodic boundary conditions in the z direction which 
are not physically appropriate. According to f|T3|) . as the film surface is approached from 
within, b± (r) — > 0. This leads to a discontinuity in across the film boundary, as can be 
seen from Eq. ( fl~5l) . and that discontinuity implies a surface current which does not exist. 

The consequences of this flaw in the Ansatz are surprisingly difficult to see — no clear sign 
of it appears in the results presented by Brandt in Ref. LL4|, many of which we reproduced 
independently. When we implemented that method the first suggestion of a problem came 
when we compared two calculations of the supercurrent which should have given the same 
results, namely j = —uQ and j = VxB = Vxb. An example is shown in Fig. [TJ for 
a system at fairly low mean induction. Note that the supercurrent calculated according 
to V x B actually circulates in the wrong direction for some values of z. A hint that the 
problem was the form of the z dependence in Eqs. (1TTT) - (1T31) . and not simply an error in 
our implementation as we first supposed, was that the disagreement become more evident 
as the maximum value of K z was increased. 

Our solution is to replace the cosine expansion for the z-dependence of q with an expan- 
sion in terms of Legendre polynomials of even order. Instead of Eqs. (JTTT)-(IT3l). take 

q(r) = &K x /^r^sinK ± • r ± P 2l (2z/d) (17) 

Ki,l 1 

b z (r) = J2b K±l cosK ± -r ± P 2l (2z/d) (18) 

K ± ,l 

h ± (r) = J2 W 1 ^ 1 sinK ± • r ± ~P 2 '< (2z/d) (19) 

K , ,Z 1 
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There is an additional benefit of the Legendre polynomial expansion for the accuracy of 
the calculations. A numerical scheme for z-integration which maintains orthogonality of 
the Legendre polynomials is appropriate for the iterative calculation of the b coefficients, 
namely, Gauss-Legendre quadrature. The abscissas for Gauss-Legendre quadrature are at 
zeros of P n (where n is larger than the highest order used in the Ansatz), and these zeros 
are more numerous near the film surfaces where the most rapid changes occur for b and q. 

We now present the full scheme for generating solutions to the GL equations for films. 
We use (• • -)u to denote the volume average over a unit cell by numerical quadrature in 
which the z abscissas are uniformly spaced, while (■ ■ -)g is the same, except it employs 
Gauss-Legendre quadrature for the z coordinate. Angle brackets without a subscript refers 
to an analytic expression for the volume average over the unit cell. Before beginning the 
iterative calculations a set of initial Ok x x z and &k ± ; coefficients must be chosen; we will 
discuss that choice following the iteration scheme. 

For the order parameter coefficients we use Brandt's iteration scheme, without modifica- 
tion, but for completeness we include it here. Defining 

g = \Vou\ 2 /4k 2 uj (20) 

the first GL equation leads to the iteration 

_ 4((w 2 - 2u + uQ 2 + g) cos K_l • rj_ cos K z z)u 
aK±Kz {5 Kzfi + mKl + Kl)/2^ + l) [Zi) 

This is always followed by an iteration to minimize F by multiplying all the ctK x x z by the 
same factor, 

a K± K z ■■= a K±Kz (u-g- uQ 2 )u/ (w 2 )u (22) 

This step was introduced by Brandt in solving the GL equations in bulk superconductors; 
if omitted, the calculations generally do not converge. 

Next comes the iteration for the &k ± 2 • Our modification of the expansions for b and q 
require corresponding changes to the iteration scheme compared to Ref. [jjl It is convenient 
to construct some auxiliary quantities such as the stray-field expansion coefficients 

& k ± = 5_>x* (23) 
i 



(compare Ref. [L4| Eqs. (10) and (21)); a quantity that arises from d(uQ 2 ) /db-x_ A 
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D K± i = (to [Q y K x - Q x K y ] sin K ± ■ v L P 2l (2z/d)) G (24) 



(compare Ref. Q Eqs. (20) and (22)); and 



S K± i = b ^' 21 ' ( 2/ ' + !) + S h ^ 1 ' 21 ( 2/ + ^ ( 25 ^ 



l'=0 l' =1+1 



which appears in 



These 
Ref. 
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d{b 2 )/db K±l = 2S K±l /(dK ± ) 2 + b K±l /(Al + 1) (26) 

ast two expressions are rather more complicated than the corresponding Eq. (19) of 
because, unlike sines and cosines, the Pn and P' 2l are not mutually orthogonal. The 
second sum in Eq. ( 125]) is finite on account of the truncation of the expansion. 
With these definitions the revised iteration scheme is 

-2S K±l - 2D K±l - 2K ± b s K± /d + c(co) b K±l 



(27) 



K 2 J{Al + 1) + c (co) 

where the constant c and the order parameter mean (u) = Y1k ± a K ± o are included to 
stabilize the iterations (compare Ref. |l4| Eq. (23)). 

The algorithm is started with an initial guess for the aK ± ^ z and &Kjj coefficients. Con- 
vergence to the physical solutions is not guaranteed, and in fact it is essential to have good 
initial values. We have used bulk solutions 15 as initial values for ok ± o an d &k ± o ; with the 
other coefficients initially zero. One repeatedly cycles through Eqs. ( 12 ip . ( I22p . and ( 12 7p until 
F has converged to an absolute tolerance of 1 x 10~ 10 or better, which typically requires 
about 200 iterations. This is slower convergence than is achieved with the cosine Ansatz 
for the z dependence for the supervelocity. A possibly related matter is that we have not 
found a suitable expression for the "mixing parameter" c that works well — large enough 
to maintain stability of the iteration scheme, small enough to allow for reasonably quick 
convergence — over the entire range of parameters that we have studied. What we do instead 
is to adjust c during the iteration cycle by monitoring the evolution of F mag and -F s tray 
We have found when either of those field energies increases excessively it is a sign that an 
instability is developing. A scheme that works reliably is that when either F mag and -F s tray 
increases by more than 50% following ( |27|) then c is multiplied by 10 and the b- iteration is 
re-run; independently, every 30 iterations c is divided by 2. 
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Although our calculations do not converge as rapidly as those reported in Ref. 14 they 
always lead to solutions with lower free energies, typically by 0.5% or less (with the same 
number of coefficients included in both calculations). These small differences are enough 
to produce noticeable changes in the phase boundaries. Our calculations also have the 
appealing feature that increasing the / cutoff for the 6k ± z always gives an improved solution; 
the same is not true of increasing the K z cutoff for the &K ± if z - Repeating the calculations 
presented in Fig. [1] yields supercurrent densities from — uQ and V x B which are nearly 
coincident, and which are close to the —uQ values displayed in that figure. 

All of the results presented in the following sections are for calculations at k — 0.5 and 
0.25; even with just those two values for the GL parameter some trends with decreasing k 
are evident. Calculations at small k are considerably more challenging: we have not yet 
been able to obtained converged solutions at k = 0.1. 



III. PHASE DIAGRAMS AND PHYSICAL PROPERTIES 

We have carried out a series of calculations at various values of B and d, and for several 
kinds of vortex lattices including triangular, square, rectangular (at various aspect ratios) 
and two classes of oblique lattices which we will refer to as rhombohedral (which interpolate 
between triangular and square at fixed unit cell area, maintaining equality of the primitive 
vector lengths) and sheared-triangular (which interpolate between triangular and rectan- 
gular at fixed unit cell area, keeping one primitive vector fixed). The common feature of 
the structures considered is that they have one vortex per unit cell, and consequently the 
coefficients in the expansion of the order parameter (fTUI) are known for a bulk system just 
below the upper critical field.— The vortex structure with lowest free energy turns out to be 
either triangular, square, or rectangular. 

Figs. |2] and [3] show the resulting phase diagrams for k = 0.5 and 0.25. The phases found 
at the upper critical field extend to lower fields, but with the phase boundaries shifting to 
larger thicknesses as B is reduced. At sufficiently low B the interval of thickness where the 
square lattice is stable is seen to vanish on the k = 0.25 phase diagram; and the same almost 
certainly holds for k = 0.5, but at a lower value of B than we have considered. Contours of 
constant aspect ratio within the rectangular phase are shown as dotted lines. On the k = 0.5 
phase diagram we have included a dashed line where we speculate that the rectangular phase 
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ends and more complicated structures with more than one flux quantum per unit cell begin; 
in drawing that line we are assuming that the aspect ratio within the rectangular phase is 
constant at the boundary with the adjacent phase. 

Within linearized GL theory dB l l 2 is constant on every phase boundary. 5 That is not a 
terrible approximation, but the numerical results are noticeably different, with the domain 
of stability of the triangular phase reduced compared to the linearized GL theory. The 
critical endpoint for the square to rectangular transition is a qualitative feature that only 
emerges from the full GL treatment. 

It is of interest to look at the free energies that underlie the phase diagrams, to see the 
scale of the free energy differences. In the lower panel of Fig. HJ F is presented as a function 
of mean induction for k = 0.5 and d = 2.4, while Fig. [5] does the same for k = 0.25 and 
d = 0.94 (the latter thickness is chosen so that the phase transitions in the two figures are 
at roughly the same values of B/k). The rhombohedral lattice free energies, not shown in 
those figures, are nearly degenerate with the free energies of square and triangular lattices at 
phase transition between them, and close to the transition their free energies almost linearly 
interpolate between square and triangular lattice free energies. 

Shear moduli have been evaluated for the three lattice structures which appear on the 
phase diagram: see the upper panels of Figs. H] and [5j For triangular lattices the only shear 
modulus is Cqq = |(cn — c 12 ). For square lattices there are two distinct types of shear, 
with moduli Cqq and |(cn — Ciz): the former preserves equality of primitive lattice vector 
length, while the latter preserves orthogonality of primitive lattice vectors. We present both 
on the figures because the latter vanishes at the continuous square-rectangular transition 
and the former is anomalously small at the discontinuous triangular-square transition. For 
the rectangular lattices we considered only the shear mode which preserves orthogonality of 
primitive lattice vectors; the corresponding modulus is |((cn + €22)/% — C12). In every case 
the shear modulus is calculated by evaluating the energy difference between the reference 
lattice structure and a slightly sheared lattice. One can see in the figures the small domains 
of metastability for the triangular and square lattice phases. It is also apparent that the 
vortex lattices at these values of k and d are anomalously soft for a wide range of mean 
inductions. 
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IV. FREE ENERGY DECOMPOSITIONS 



The preceding section presented the main physical results of the calculations, but further 
insight might be gained by comparing not just F for different lattice structures but also 
various "components" of the free energy. 

One decomposition is into the condensation, kinetic, and magnetic terms described fol- 
lowing Eq. ([3]). Let us first consider k = 0.5, B/k = 0.825, and d = 2.0, which is in the 
triangular phase but not far from the square phase. For a bulk system at the same GL 
parameter and mean induction, the square lattice has lower free energy density than the tri- 
angular lattice. Why is the relative stability reversed? In Table [J we present the differences 
in free energy density components between the film and the bulk system for both triangular 
and square vortex lattices. The signs of all those differences may be understood as a con- 
sequence of suppression of the order parameter in the film compared to the bulk. However, 
the exchange of stability is a more subtle matter, since that depends on the difference (tri- 
angular minus square lattice values) of those free energy density differences. Alternatively, 
we can compare the triangular and square lattice free energy density components for films of 
different thickness, as presented in Table HT1 It is then evident that with increasing thickness, 
the transition to the square vortex lattice is favored only by the condensation term. 

We can also examine the z-dependence of the free energy density, integrating in Eq. ([3]) 
only over x and y and dividing only by S to define F(z). (F stray is taken as a ^-independent 
contribution to F(z).) Figure [6] compares square and triangular vortex lattices for k = 0.5 
just below the upper critical field for d = 1.5, 2.0, and 2.5. The triangular lattice has lower 
total free energy only for d — 1.5. However, in every case F(z) is lower for the triangular 
lattice when z ~ d/2, and, with decreasing z, F(z) decreases more rapidly for the square 
lattice than for the triangular lattice. Figure [6] is thus consistent with the interior of the 
film being more bulk-like than the surface; and in fact F(0) approaches F for a bulk system 
with as d increases. In terms of the free energy components, the condensation term is nearly 
independent of z, as is uj (as pointed out by Brandt for films of type-II superconductors^). 
The kinetic term is responsible for the z-dependence seen in Figure O since the magnetic 
term is smaller at the surface, where the field lines spread out, than in the center of the film. 
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V. CONCLUSIONS 



We have improved Brandt's methocU 4 - for solving the GL equations for thin-film super- 
conductors in perpendicular magnetic fields, and applied it to a series of calculations for 
various vortex lattice structures with one vortex per primitive cell in type-I superconductor 
films of intermediate thickness {d ~ A). The phase diagrams presented in Sec. (IHip are 
the first step beyond the linearized theory towards the development of an accurate equi- 
librium flux structure phase diagram for films of type-I GL superconductors. The results 
suggest that non-triangular flux lattice structures (square and rectangular) may arise at 
mean inductions well below the upper critical value. For future work we wish to carry out 
similar calculations for flux structures with more that one vortex per primitive cell, as well 
as lattices of multiply-quantized vortices and various intermediate state models; these will 
require different expansions for the in-plane variation of u, q and b but Legendre function 
expansions for the z dependence should still be applicable. 

The anomalous softness of the vortex lattice in and near the domain of stability for 
the square vortex lattice offers hope that some features of the theoretical phase diagram 
might be observed in critical current measurements, in the form of a "peak effect"— well 
below the upper critical field. However, quantitative comparison between the theoretical 
phase diagrams and experimental results will necessarily be complicated by anisotropy and, 
possibly, thermal fluctuations.— 
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Table I: Film minus bulk free energy density terms for k = 0.5, B/k = 0.825, and d = 2.0. 





Triangular 


Square 


10 4 AF cond 


161. 


180. 


10 4 AF kin 


-116. 


-129. 


10 4 (AF mag + AF stray ) 


-25.5 


-29.0 



Table II: Differences in free energy density terms (triangular lattice minus square lattice) at k = 0.5 
and B/k = 0.825 for several values of d. Note that AF < for d = 2.0 but is positive at the other 
thicknesses. 



d 


2.0 


2.33 


2.6 


10 4 AF cond 


-4.44 


-2.64 


-1.37 


10 4 AF kin 


6.55 


5.29 


4.42 


10 4 (AF mag + AF stray ) 


-2.31 


-2.60 


-2.81 
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Figure 1: Supercurrent density components j x and j y calculated from — cjQ (crosses) and V x B 
(circles) from solutions by the method of Ref. [l4j for a system with k = 0.5, B = 0.4/k, d = 4.3, 
and a 32 x 13 x 9 grid for real-space sampling. The vortex lattice is triangular, with one primitive 
translation being x\x. In these plots y = 0.017xi, with z = 0.89d/2 for (a) and (b) and z = for 
(c) and (d). 
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Figure 2: Vortex lattice phase diagram for k = 0.5. The triangular-square transition is discontinuous 
while the square-rectangular transition is continuous. Inside the rectangular phase, the dotted lines 
labeled 0.6 and 0.4 are contours of constant aspect ratio. The dashed line corresponds to the aspect 
ratio of 0.36, which is the smallest aspect ratio for which a rectangular lattice is stable at the upper 
critical field (B/k = 1), following Callaway.— 
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Figure 4: Shear moduli and free energies per unit volume for triangular, square, and rectangular 
lattices, at k = 0.5 and d = 2.4, for mean inductions around the domain of stability of the square 
lattice. Free energies are referenced to value for the square lattice; on that graph the triangular 
lattice values are the triangles and the minimum- F rectangular lattice values are the circles. The 
vertical dashed segments in both plots indicate the transition between triangular and square lattices, 
to make clear the discontinuity in shear modulus. On the shear modulus plot, triangles are cqq for 
the triangular lattice, diamonds are cqq for the square lattice, squares are \ {c\\ — cu) for the square 
lattice, and circles are \{(c\i + 022)/^ — C12) for the minimum-F rectangular lattice. Both F and c 
are in units of HqH^. 
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10 5 F(z) 

Figure 6: Free energy density dependence on z for films of various thickness with k = 0.5 and 
B = 0.99/re, for square and triangular vortex lattices (indicated by the squares and triangles, lines 
are guides for the eye). 
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